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Abstract 

We determine the nuclear modifications of parton distribution functions of bound pro- 
tons at scales Q 2 > 1.69 GeV 2 and momentum fractions 10~ 5 < x < 1 in a global 
analysis which utilizes nuclear hard process data, sum rules and leading-order DGLAP 
scale evolution. The main improvements over our earlier work EKS98 are the au- 
tomated x 2 minimization, simplified and better controllable fit functions, and most 
importantly, the possibility for error estimates. The resulting 16-parameter fit to the 
N = 514 datapoints is good, xVd.o.f = 0.82. Within the error estimates obtained, the 
old EKS98 parametrization is found to be fully consistent with the present analysis, 
with no essential difference in terms of x 2 either. We also determine separate uncer- 
tainty bands for the nuclear gluon and sea quark modifications in the large- 2 region 
where they are not stringently constrained by the available data. Comparison with 
other global analyses is shown and uncertainties demonstrated. Finally, we show that 
RHIC-BRAHMS data for inclusive hadron production in d+Au collisions lend support 
for a stronger gluon shadowing at x < 0.01 and also that fairly large changes in the 
gluon modifications do not rapidly deteriorate the goodness of the overall fits, as long 
as the initial gluon modifications in the region x ~ 0.02 — 0.04 remain small. 
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1 Introduction 



Universal, process-independent parton distribution functions (PDFs) of free and bound 
nucleons are a key element in the computational phenomenology of processes involving 
large virtualities Q 2 in hadronic and nuclear collisions. The free proton PDFs are 
nowadays rather well constrained through the global analyses [H El E] , which use the 
DGLAP jl] Q 2 -evolution, sum rules and a large amount of data from deep inelastic 
lepton-proton scattering (DIS) and high energy proton-(anti)proton collisions. The 
success of the forthcoming Large Hadron Collider (LHC) program in the search for the 
Higgs boson and physics beyond the Standard Model depends on the precision of the 
PDFs. 

At collider energies, hard processes are abundantly available also in heavy-ion colli- 
sions. These processes play an important role in testing QCD dynamics and factoriza- 
tion, as well as in the search of quark-gluon plasma signatures and in the determination 
of the QCD matter properties. Similar to the free proton case, the computation of nu- 
clear hard process cross sections requires the nuclear parton distributions (nPDFs) as 
input. Thus, there is an obvious need for the global analyses of the nPDFs such as 
presented in [51 El El El E] ■ 

Hard partonic processes taking place at mid-rapidities in nuclear collisions at the 
Relativistic Heavy Ion Collider (RHIC; A+A and d+Au at -\/s NN = 200 GeV) typically 
probe the nPDFs in a kinematic region where the nuclear effects remain relatively small 
and are fairly well constrained by the global analyses. Towards smaller scales and off 
mid-rapidity, however, the probed region extends towards smaller momentum fractions 
x where both the nuclear effects and the uncertainties in the nPDFs grow larger. Soon 
at the Large Hadron Collider (LHC; Pb+Pb at \fs NN = 5.5 TeV) the range of scales 
and fractional momenta probed will be widened further, both towards smaller x and 
towards larger Q 2 . This, together with the fact that the nuclear gluon distributions 
are still relatively badly known, emphasizes the importance and topicality of the global 
analyses in pinning down the nPDFs and their uncertainties. 

The fact that nuclear and free proton PDFs are mutually different has been known 
for well over twenty years; for a recent review, see Ref. [10J. The nuclear effects, the nu- 
clear modifications relative to the free proton PDFs, are usually named according to the 
observed behaviour of the nucleus-to-Deuterium ratio of the structure functions in 
different x-regions, as follows: (i) shadowing; a depletion at x<0.1, (ii) antishadowing; 
an excess at 0.1^x^0.3, (iii) EMC effect; a depletion at 0.3<x<0.7 and (iv) Fermi 
motion; an excess towards x — > 1 and beyond. This nomenclature will be used in this 
paper as well. The dynamical origin of these nuclear modifications has been actively 
studied in different frameworks as well, see the Refs. e.g. in [TOl HH U2] . The DGLAP 
evolution of the nPDFs and their modifications relative to the free proton PDFs have 
been studied for two decades, see e.g. Refs. HH [151 HH H71 HH UHl [201 [2U H2] . 

In a global DGLAP analysis the nPDFs are pinned down as model-independently 
as possible at a chosen initial scale on the basis of DGLAP evolution, sum rules and 
hard process data from nuclear collisions. So far, three groups have presented global 
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DGLAP analyses of the nPDFs analogous to those of the free proton. These are the 
ones by us, Eskola et al. EKS98 [HIE], by Hirai et al. HKM [7J and HKN [8], and by 
de Florian and Sassot nDS [H] . The EKS98 analysis [5J E] was the first one to show 
that a good overall fit to the nuclear DIS and Drell-Yan (DY) data can be obtained 
in a DGLAP-based global analysis. In particular, the scale-dependence of the ratio 
F2/F2 observed by the NMC experiment jT3j was very nicely reproduced by tuning 
the initial gluon modifications suitably. The iterative x 2 minimization in EKS98 was 
carried out manually (by eye), and no well-controlled error estimates were obtained. 
Since then, extensive further work has been done by Kumano and his collaborators in 
estimating these uncertainties [7J [8], and by de Florian and Sassot [9] in bringing the 
global nPDF analysis to the next-to-leading order (NLO) level. 

In this paper, we perform a global analysis of the nPDFs in the EKS98 framework. 
Our study is partly a reanalysis of EKS98 as we take some guidelines from this old fit. 
To minimize the number of fit parameters, however, we now apply simpler piecewise 
analytical shapes for the nuclear effects at the initial scale. We also construct the 
nuclear quark modifications in a more transparent way than in our previous work. 
The goal here is twofold: on one hand, by making the x 2 minimization procedure 
automated, we wish to check whether the goodness of the old EKS98 fit could still be 
improved, and on the other hand we wish to get a better hold on the uncertainties of 
the nPDFs, of the gluons in particular, in this framework. 

The results of this study can be summarized as follows: Within the obtained x 2 an d 
error estimates, we conclude that the old EKS98 parametrization still serves very well. 
Thus, we do not release a new parametrization but recommend to use EKS98. We 
also demonstrate how the small-x nuclear gluon distributions are, in spite of the good 
overall fit obtained, still not well constrained with the currently available nuclear DIS 
and DY data elsewhere than perhaps at x ~ 0.02 — 0.04. A comparison with the results 
from the previous global analyses is also shown, demonstrating the nPDFs uncertainties 
concretely. Finally, a special case beyond the original EKS98 setup, a gluon shadowing 
clearly stronger than that in F^/F®, is considered and further developments of the 
analysis by inclusion of RHIC data are discussed. 

This paper is organized as follows. In Sec. [21 we define nPDFs according to the 
EKS98 framework and introduce the fitting procedure. Section [3] contains the results 
of x 2 minimization and the a detailed comparison with the nuclear DIS and DY data. 
Section H] is devoted for the comparison with previous global analyses. In Sec. [5] we 
show the results from the error analysis performed and verify the validity of EKS98. 
In Sec. El we discuss the possibility of a stronger gluon shadowing supported by the 
RHIC data. Conclusions and further discussion are given in Sec. [71 
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2 The framework 



2.1 Definition of nPDFs 

As introduced in EKS98 |5], by a nuclear parton distribution function f A we refer to 
the distribution of a parton type % in a proton^ bound to a nucleus of a mass number A. 
We define and parametrize the nuclear modifications relative to the known free proton 
PDFs 

In the EKS98 framework which we adopt here, the PDFs of the bound neutrons are 
obtained from f A (x,Q 2 ) by assuming isospin symmetry. Thus, e.g. the total w-quark 
distribution in a nucleus of a mass number A and a proton number Z becomes Ua = 
Z fu + {A ~ Z)fd- Correspondingly, the lowest-order QCD parton model expression 
for the I A DIS structure function F 2 then becomes F A = J2q 6 q [Qa + Qa] > where 
Q = U,D,S,.... 

The total amount of fit parameters in the initial ratios R A must be limited for ob- 
taining converging well-constrained fits. Unfortunately, the variety of the nuclear data 
is presently not enough to pin down each R a (x,Qq) separately. Therefore, following 
the EKS98 procedure, we can include only three different ratios for each nucleus at an 
initial scale Q 2 = Qq where heavy quarks can be neglected: The same average modifi- 
cation Ry = (f A v + fd v )/(fu v + fd v ) is applied for all valence quarks separately (only 
at Ql however), the corresponding sea quark average modification R A applied for all 
sea quarks separately (again at Ql only) and Rq for gluons. While this is the best we 
can do here, we note that the valence u and d quark nuclear modifications may in fact 
well differ from each other - for a recent study of how large differences between R A 
and R A would explain the NuTeV weak-mixing angle anomaly observed in v{v)+~Fe 
DIS, see [23J. Also, in the sea quark sector, due to their mutually differing absolute 
distributions, it would be natural to expect that the initial s quark modifications are 
not necessarily identical to those of u and d. Without a multitude of further data 
constraints, however, such details cannot be reliably included in a global analysis. 

In the original EKS98 analysis [5] we first parametrized the DIS structure function 
ratio 

at the initial scale and then decomposed this into the valence and sea parts. The 
initial gluon modifications were obtained by adding a double gaussian distribution on 
the antishadowing peak of the parametrized Rp . In the current analysis we choose a 
more straightforward procedure by parametrizing directly the ratios Ry, R A and Rq 
at Qq- 

1 Note that in HKN a slightly different definition is used, see [7J [5] . 
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The initial scale is here chosen to be Qo = 1.3 GeV in order to match the CTEQ6L1 
PDF set [3], which we use to calculate the absolute nuclear PDFs at Qq\ 

ff(x,Ql)=R?(x,Ql)fr m6L1 (x,Ql). (3) 

The lowest order DGLAP scale evolution is calculated using the routine from the CTEQ 
collaboration [22] as it provided fast enough evolution for the minimization purposes. 

The key constraints for the nPDFs are given by the nuclear hard process data from 
lepton-nucleus DIS and from the DY dilepton production in proton-nucleus collisions. 
We utilize the results from the DIS measurements, available in the form of ratios over 
Deuterium and Carbon, 

\da lA /dQ 2 dx l _q a \da lA /dQ 2 dx L p Rf. 2 (x, Q 2 ) 

\da m /dQ 2 dx~ h ±da lc /dQ 2 dx~ R%{x,Q 2 Y U 

where the LO connection is implied. The DY data are available in the form of ratios 
over Deuterium and Beryllium, 

\do- p D y/dx 2 dQ 2 lp . gjggy / 'dxxdQ 2 l _q R% y {x u Q 2 ) 

\daf Y /dx 2 dQ 2 Kd ^ X2 ^>> \daf«/d Xl dQ 2 Rfy^Q 2 )' W 

Above, Q 2 is the invariant mass of the dilepton pair and Q 2 = X\X2\fs NN . The data 
included in this study are shown in Table [H The small nuclear effects in Deuterium 
are neglected. 

As will become clear in the error analysis presented in Sec. [5j the available sets of 
experimental data do not constrain the distributions of different parton flavours over 
the whole range of x. This will be reflected as some assumptions regarding the shape of 
the ratios which are basically the same as in our previous EKS98 work. In particular, 
motivated by the requirement of a stable evolution (that the nuclear modifications 
should not change very rapidly from their starting values), a saturation (flattening) 
of the ratios R A at x — > 0, and a valence quark -like behavior of the sea and gluon 
modifications for x — > 1 will be assumed. In the following we explain in detail how the 
initial parametrization for R A (x, Qq), R a (x, Qq) and R A (x, Qq) was constructed. 

2.2 Fit functions and parameters 

While the basic idea in the global DGLAP analysis is straightforward, it is a surpris- 
ingly nontrivial task to develop functional forms for the fit functions for the ratios 
Ry, R A and Rq which can be used in the automated \ 2 minimization process in a 
transparent way. To have a better control over the multidimensional parameter space 
and over the numerical results obtained, each parameter should preferably have a clear 
interpretation, too. Due to the various A and x dependent nuclear effects discussed 
above and also due to the mutual differences between the valence, sea and gluon modi- 
fications, the fit functions must contain sufficiently many parameters to secure enough 
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Table 1: The data used in this analysis, grouped according to the nuclei measured. The 
mass numbers are given in parentheses. The number of datapoints refers to those falling into 
the region Q 2 > Q%. 



5 



flexibility necessary for obtaining good fits. At the same time, the number of parame- 
ters has to be reduced to a minimum in order to obtain converging fits with the rather 
limited set of data constraints at our disposal. Finally, once the working functional 
forms have been verified, one needs to analyze (on the basis of the data constraints 
and x 2 fits) which parameters can be left free and which can be fixed. Furthermore, 
the best local minimum in \ 2 has to be verified by optimizing the the initial values of 
all free parameters. All this implies extensive manual labour, even though the actual 
search for the x 2 minimum is automated. 

For the controllability discussed above, and after various other attempts, we ended 
up constructing each of the initial ratios R v , R A and Rq from three different pieces: 
R A (x) at small values of x below the antishadowing^] maximum, x < x A ; R A (x) in the 
medium- x region from the antishadowing maximum to the EMC minimum, x A < x < 
x A ; and R A (x) in the Fermi-motion region in the large- x region, x > x A - 



R A (x) = c a + (c A + c A x)[exp(~x/xf) - exp(-xf /x A )}, x < x A (6) 
R A (x) = a A + a A x + a A x 2 + a A x 3 , x A < x < x A (7) 



x) 



3 w ft ~,\/3 A ' e 



In choosing the above forms, we were motivated by the functional forms used before 
in Hard Probes [21] (see [32]), EKS98 [5] and HKN [S]. Matching is done by requiring 
continuity of the fit functions and setting their first derivatives to zero at x A (local 
maximum) and x A (local minimum). As the coefficients af , bf and cf are somewhat 
unintuitive, we shall quote the results in terms of the following more transparent set 
of seven parameters from which these coefficients can be easily solved: 
y£ Rf at x — > 0, 

slope factor in the exponential, 
x a ? Va position and height of the antishadowing maximum 

He position and height of the EMC mimimum 
j3 A slope of the divergence of R3 at x — > 1. 

Each of the above parameters is in principle yet specific to a nucleus A. This (at 
least) doubles the amount of parameters. We parametrize the A-dependence in a simple 
power-like form: 

^ = ^ rcf (^-) Pz S (9) 

where Zi = x s , x a , y a . . ., and choose the reference nucleus to be Carbon, A re { = 12. 
The number of parameters we have for the valence, sea and gluon ratios each is thus 
14: the Carbon parameters (suppressing the superscript C to lighten the notation) y Q , 
x s , x a , x e , y a , y e , p, and their powers p yo , p Xs , p Xa , p Xe , p Va , p Ve and pp. Altogether 
this makes 3 x 14 = 42 free parameters. Even if the momentum and baryon number 



2 For antiquarks Rg < 1, by antishadowing we refer to the shape similar to Rp 2 - 
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conservation, imposed individually for each nucleus, reduce this number by four, it is 
clearly far too large for a converging \ 2 minimization process, given the limited data 
constraints we have. In order to radically reduce the number of free parameters, we 
proceed as follows, keeping in mind the focus on the small- and medium- a; regions. 

• Fermi-motion. In the large- x region, where valence quarks dominate, the DIS or 
DY data do not give proper constraints for gluons or sea quarks. Thus, we fix the 
Fermi-motion slopes (3 A in R A and Rq to be the same as in R A . Based on our 
previous EKS98 work, we fix (3 = 0.3 and pp = in R v , thus ignoring a possible 
A-dependence of (3 A . 

• EMC effect. Gluons originate from valence quarks at small scales and large x. 
Therefore, they should reflect the EMC effect observed in Ry (Rp ). From the 
gluons the effect should then be transmitted on to R A as well. We have checked 
that this is indeed the case in the DGLAP evolution [33]. Thus, by assuming the 
similarity of the EMC-minima in each initial ratio Rq, R a and R A , one reaches 
a stable scale evolution of this nuclear effect. As the available data, however, 
constrain the EMC effect in detail only in R A , we fix the location parameters 
x A and the magnitude parameters y A of the EMC-minima in Rq and R A to be 
identical to those in Ry. For the valence part, we noticed that allowing for an 
A dependence in x A did not improve the overall fits, hence we fix p Xe = for 
simplicity. 

• Antishadowing. In course of the present analysis we also noticed that the location 
parameters x A of the antishadowing maxima in R A and in R A typically become 
almost A-independent and that the weak A dependence does not improve the 
obtained fits. We therefore set p Xa = in R A and R A . In order to reduce the 
number of gluon parameters to the very minimum, we simply fix x A of gluons to 
be identical to that in valence but leave y a and p yA free for controlling the height 
of the antishadowing maximum in an A-dependent way. 

• Shadowing. In the small-a: parts, based on x 2 -checks, we drop the A-dependence 
of the slope parameters x A , hence setting p Xs to zero and and leaving x s free in 
all ratios. 

• Conservation laws. Baryon number and momentum conservation are used to 
calculate y A for R A and Rq, respectively, for each nucleus individually. This 
eliminates the parameters yo and p yo for the valence and gluon modifications. 
For the sea quarks, these parameters are left free. 

All this brings the number of free parameters down to 16: x s , x a , x e , y a , p Va , y e 
and p ye in R A ; y , p yo , x s , x a , y a and p ya in R A ; and x s , y a , and p Va in Rq. Table [2] 
summarizes the above discussion on the parameters as well as their values obtained in 
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finding a "best" local minimum with respect to the fit parameters for 



A^data / 



dataj — theory ^ 



) 



2 



i=l V 



A, 



(10) 



As the data errors Aj, we take the given statistical and systematic errors added in 
quadrature. 

Some remarks on the functional form adopted for the shadowings at small- £ are in 
order here. Since the valence modification R v is rather well constrained by the DIS and 
DY data in the large- and medium-x regions, its small-x behaviour becomes relatively 
stringently constrained by the baryon number sum rule. Unfortunately, in the absence 
of DIS (or DY) data for R F2 at x < 0.001 in the DGLAP region Q 2 > 1 GeV 2 , the sea 
quark R$ and the gluon Rq cannot be pinned down similarly well in the small- x region 
- thus their behaviour and error estimates at small x are bound to be specific to the 
fit function forms assumed. 

The motivation for choosing the smallest-x form of Rf(x) in Eq. where shadow- 
ing levels off to a constant value at x = 0, is the fact that such saturation of shadowing 
has been observed in the very small-x & very small-Q 2 DIS data (see Fig. 10 in [28J) 
and the fact that the Q 2 dependence there is rather weak (see Figs. 11 and 12 in |28j). 
In doing this, however, we should keep in mind that the implications of the observed 
saturation of shadowing are not clear for the nPDFs at perturbative scales: power 
corrections ~ {Q 2 )~ n [31] are most likely important in the DIS cross sections at small 
enough scales, and also nonlinearities [35] (neglected here) are expected to play a role 
in the scale evolution at sufficiently small-x & small-Q 2 . 

In the previous EKS98 analysis, due to the modest and non-negative log Q 2 -slopes 
of Rp 2 discussed above, we fixed the smallest-x behaviour of Rp 2 (x, Ql) to a value 
slightly above the saturation of shadowing observed at lower scales. The log Q 2 slopes 
of Rp 2 computed from the DGLAP equations at small x [311 El E2], 



are non-neg ative if R&(2x) > Rp 2 (x). In EKS98, it was shown that an ansatz Rq(x — > 
0) — > Rp 2 (x — > 0) works well for the smallest x. In the present analysis, we want to test 
the above EKS98 gluon framework and thus keep the saturation of gluon shadowing 
independent of that in Rg. 

3 Results 

3.1 Final parameters and their interpretation 

In minimizing the x 2 with respect to the free parameters, we used the MINUIT routines 
from the CERN Program Library [37] . Only after reducing the number of free param- 
eters down to 16, and after extensive searches for suitable initial parameter values, we 



dR F2 (x,Q 2 ) 
d log Q 2 




(11) 



8 



were able to find a converging fit indicating a local minimum of the x 2 ■ The obtained 
parameters for the best fit found are shown in Table El The resulting goodness of the fit 
was x 2 = 410.15 for N = 514 data points and 16 free parameters, giving x 2 /N = 0-80 
and xVd-o.f. = 0.82. 
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Table 2: List of all parameters defining the modifications Ry, Rg and Rq in Eqs. ((MHl) 
at the initial scale Qq = 1.69 GeV 2 . The parameters yo, y a , He, x a , x e and j3 are for 
the reference nucleus A = 12, and the powers pt define the ^4-dependence in the form of 
Eq. 0. The obtained final results for the fitted 16 free parameters are shown and the fixed 
parameters are indicated. The parameters which drifted to their upper (u) and lower (1) 
limits are indicated, see the text for details. 

As indicated in the table, the parameters x s controlling the slopes of R± near the 
antishadowing region were drifting to their limits. In spite of various attempts we failed 
to improve upon this unwanted feature. Obviously, there is still room for developing 
the chosen functional forms in the quark sector too. However, as the fits obtained 
now (and already in EKS98) are very good, new functional forms are not likely to 
improve the x 2 essentially. In fact, this was our observation also at different stages of 
the present analysis: in spite of the non-converging fits often obtained (which were due 
to too many free parameters allowed or badly guessed initial parameter values), the 
obtained fits themselves were equally good. 

The gluon sector, however, is the most troublesome one, as all the data constraints 
are indirect and not very conclusive when put into the context of a global analysis: 
rather large changes in the gluon shadowing and antishadowing can be compensated 
for by fairly moderate modifications in the quark sector. As a result, gluons have a 
minor effect in the overall x 2 ■ The gluonic parameters y a and p ya , which are drifting 
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to their limits (see the table), reflect these problems. 

As described above, the functional form R± at very small x preassumes the satu- 
ration of shadowing also for gluons. The height of the antishadowing bump y a and its 
A-dependence are correlated with the parameters yo and its A dependence p Va which 
are computed from the momentum sum rule: the larger the y a , the smaller the yo. 
Even though no essential improvement over the \ 2 was noticed in varying the limits of 
y a and p Va , a clear trend was observed: as indicated by reaching the lower limit of y a , 
the amounts of gluon antishadowing and shadowing always tend to be minimized. This 
in turn means that gluon shadowing saturates at a value larger than that of sea quarks 
and that the log Q 2 slopes of Rp 2 at the smallest x remain positive. These observations 
coincide with the results from previous global analyses HKN [8] and nDS [9] . 

We thus conclude that the present DIS and DY data and the sum rule constraints 
suggest that gluon shadowing is weaker or at most as strong as that in sea quarks. As 
one of the goals here is to test the EKS98 framework for our final results summarized 
in Table [2] we have set the lower limits of the free gluonic antishadowing parameters y a 
and p Va in such a way that the gluon shadowing levels off to the same value as that of 
sea quarks (Rp 2 ). The benefit in doing this is that we can keep the EKS98-\ike good 
agreement with the clearly positive logQ 2 slopes of F 2 Sn /F 2 c observed at x ~ 0.01, see 
Fig. [9] ahead. 

As explained above, in the present analysis the valence and gluon parameters y$ are 
computed from baryon number and momentum sum rules, correspondingly, for each 
nucleus separately. For completeness, we note that a power-law fit of Eq. ([9]) to the 
values obtained, using A = 12 and 208, gives y = 0.9288 and p yo = —0.031209 for 
valence and yo = 0.8898 and p yo = —0.084315 for gluons. With such parametrization, 
baryon number and momentum would be conserved with sufficient accuracy, within a 
few per cent, for all nuclei. 

The obtained initial nuclear modifications are shown in Fig. [H where we plot 
Rv( x i Qo) ( s °lid lines), R$(x, Ql) (dotted lines), Rq(x, Ql) (dashed lines) and Rp 2 {x, Ql) 
(dotted-dashed lines) for nuclei A = 12, 40, 117 and 208 at an initial scale Ql = 1.69 
GeV 2 . 

The scale evolution of the nuclear effects is shown in Fig. [2j where the ratios are 
plotted for A = 12 and A = 208 as a function of x, at fixed scales Q 2 = Ql = 1-69 GeV 2 
(solid), 10 GeV 2 (dotted) and 10 4 GeV 2 (dashed). In the regions where no stringent 
data constraints are available for sea quarks and gluons, notice the systematic scale 
dependence at small x (logQ 2 slopes do not change their sign), and the stability of the 
ratios near the EMC minimum. 

3.2 Comparison with data 

Next we compare the obtained results with the data included in the analysis and 
illustrate the good overall agreement obtained. The DIS data can be found in Figs. [3H6] 
and in[9l and the DY data in Figs. [THE! In the plots below, the statistical and systematic 
errors of the data have been added in quadrature. 



10 



1.3 
1.2 
1.1 
1.0 
0.9 
0.8 
0.7 

o 

Q> 

S 1.3 
1.2 
1.1 
1.0 
0.9 
0.8 
0.7 
0.6 



A=12 












1 1 1 Mill 


— 1 1 1 Mill 


— 1 1 Mill 


— 1 



R F 2 

I I I M il 



A=m 



~ rTTTT ] rn 

A=40 









1 1 1 Mill 


— 1 1 1 Mill 


— 1 1 1 Mill 



I I M ill 



,4=208 




1.3 
1.2 
1.1 
1.0 
0.9 
0.8 
0.7 

1.3 
1.2 
1.1 
1.0 
0.9 
0.8 
0.7 
0.6 



10" 



10 



-3 



10" 



10" 



10" 



10" 



10" 



10" 



1 



X 



x 



Figure 1: Initial nuclear ratios Ry(x, Qq) (solid lines), Rg(x, Qq) (dotted lines), Rq(x,Qq) 
(dashed lines) and R§, (x,Q%) (dotted-dashed lines) for A = 12, 40, 117 and 208 at Ql = 1.69 
GeV 2 . 

In Fig. [3] we show the computed ratio -^FA/-Lf£ = Rp 2 / R% against the NMC 
data [29J for various nuclei. The open squares are the NMC data points and the filled 
squares are our results computed at the corresponding values of x and Q 2 . This data 
set plays a major role in constraining the A-systematics of nuclear quark distributions 
at small x. 

In Fig. H]we compare the computed ratio Rp 2 (x, Q 2 ) with the data from SLAC [25] . 
E665 [26], NMC 95 [28] and NMC 95 [27] reanalysis. The open triangles, diamonds, 
squares and circles stand for the data and the corresponding filled symbols show our 
results. Note that at the same/similar values of x the values of Q 2 can vary between 
the different data sets, hence the multiple filled symbols at these x. In the figure, we 
have also included the small-x data points whose Q 2 -values lie below our initial scale. 
The asterisks show our results at our Qq. To compare these points with the data, 
one should perform the scale evolution downwards. We do not consider this here (and 
hence these data points are not included in the \ 2 minimization either) but from the 
figure we can immediately see, as the logQ 2 slopes of Rp are positive and modest, and 
as the points computed at a higher scale lie above the NMC data, that the agreement 
is good also in that part of the small-Q 2 region where the DGLAP might still be valid. 
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Figure 2: Scale evolution of nuclear modifications: the ratios R v (x,Q 2 ), R v (x,Q 2 ), 
R v (x,Q 2 ), and Rj±(x,Q 2 ) at scales Q 2 = 1.69, 100 and 10000 GeV 2 for A = 12 and 208. 
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Figure 3: The computed ratio Rp 2 (x, Q 2 ) vs. Rp 2 (x, Q 2 ) compared with the NMC data [29] . 
The open symbols are the data points with errors added in quadrature, the filled ones are 
the corresponding results from this analysis. 
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Figure 4: Calculated Rp 2 (x,Q 2 ) (filled symbols) are compared to SLAC (triangles) 
E665 (diamonds) [26J, NMC 95 (squares) [25] and reanalysed NMC 95 (circles) data [27]. 
The asterisks denote our results calculated at the initial scale Qq, these are for the smallest-x 
data points whose scales lie in the region Q 2 < Qq. 
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Similar comparisons are shown in Fig. [5] for the ratios Rp^/Rp 1 and Rp 2 /Rp 2 . In 
the upper panel we show the ratio Rp? / Rp 2 from the E665 experiment (open triangles) 
|2bj . The agreement is not very good, which is not surprising as the NMC and E665 
data sets in Fig. H] do not agree, either (the NMC data has more weight in the analysis 
due to their smaller error bars). However, as noticed by the NMC well in the past [29], 
if one considers the ratio of ratios, Rpl° / Rp 2 , the agreement between these data sets 
becomes very good. This is shown in the lower panel of Fig. [5j where we plot the 
data from NMC (open squares) [29j and together with a ratio calculated from the E665 
(open triangles) data for R F ^ and Rp [26] . We obtain the error bars for the computed 
E665 Pb/C ratio by first adding the statistical and systematic errors in quadrature 
separately for Pb/D and C/D, and then taking these errors to be independent. The 
filled squares and triangles again show our DGLAP results corresponding to the data 
points, while the asterisks mark our results at the x-points where our initial scale is 
higher than the Q 2 in the E665 data. 

Further comparison with the SLAC data [25] for R$ 2 {x, Q 2 ) are shown in Fig. [6] for 
various nuclei and Q 2 scales. This set of data plays an important role in constraining 
x- and A-dependence of the valence quark distributions in the EMC region. The filled 
symbols again stand for our results, the open ones for the data. 

Figure [7| shows the comparison of the calculated LO Drell-Yan cross section ratios, 
Eq. (EJ), to the FNAL E772 data |24j . The momentum fraction x 2 is that of the nuclear 
parton. Open squares with error bars present the data points and filled squares show 
the calculated values. As can be seen, the calculated values fit the data rather well, 
except at the smallest appoints for tungsten (for which the EKS98 seems to work 
slightly better). 

Figure [8] then shows the comparison with a newer E866 data set [30] on the DY 
ratio (da pA / dQ 2 dxi) / (da pD / dQ 2 dxi) as a function of the projectile-parton momentum 
fraction. Four different invariant mass bins are considered. Large values of X\ now 
correspond to small values of x%. Confirming the trend seen in the previous figure, we 
note that the A dependence of shadowing could be slightly stronger in order to better 
match with the DY data. Within the present global analysis, however, we were unable 
to improve on this feature. 

Finally, in Fig. [9] we plot the scale evolution of the ratio -^F^ / 12^2 compared 
with the data from NMC [13] for several fixed values of x. The logQ 2 slopes of the 
data at small x, which are sensitive to the gluon modifications as shown in Eq. (TTTI) are 
reproduced very well, similar to EKS98. Note that the 15 panels here correspond to the 
15 data points in the lower left panel of Fig. [3], so that the normalization of -jj^F^* 1 / j2^2 
at each x is given by the overall fit. Thus in the upper left panel (x = 0.0125) of Fig. [9] 
the normalization is slightly higher than that of the data, while in the third panel 
(x = 0.025) both the normalization and the logQ 2 slopes match perfectly. 

The NMC data at the smallest-x panels of Fig. [9] play an important role in con- 
straining the nuclear gluon modifications. These data were the key ingredient in the 
EKS98 analysis in pinning down the nuclear gluon modifications around x ~ 0.03, for 
more discussion see also [32J. We note, however, that in an automated global analysis 
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Figure 5: Top: The ratios Rjp'/Rjp from the E665 experiment (open triangles) [26J com- 
pared with the results from the present analysis (filled triangles). Bottom: Comparison of 
the ratios R^/R%. The NMC data [29] are shown by open squares, the ratios calculated 
from the E665 data [26] by open triangles. For the error estimates in the latter case, see the 
text. The corresponding theoretical results are again shown by the filled symbols, and by 
asterisks if the experimental Q 2 is below our initial scale Qq. 
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Figure 7: The ratio of the computed LO differential Drell-Yan cross sections (open squares), 
(da pA /dQ 2 dx 2 )/{da pB /dQ 2 dx 2 ), and the E772 data [24J (filled squares). 



like we perform here, this role becomes not quite as clear: even relatively large varia- 
tions of the gluon modifications induce changes practically only in the first few panels 
of this figure. The weight that these these panels have in the x 2 is rather small among 
the 500 other data points from cross sections mostly sensitive to the changes in the 
quark sector. 



4 Comparison with previous analyses 

Table [3] summarizes the \ 2 obtained in this work, EKS98 \6\ , HKM [7], HKN j8] and 
nDS |9] analyses. Since each analysis uses different initial scales, different amount of 
data points and different data sets, we quote the values given in the original references 
(except for EKS98 whose x 2 we compute here using CTEQ6L1). As seen in the table, 
the goodness of the fit using the EKS98 nuclear effects is very close to the one obtained 
in this work and also (contrary to the claim in [9]) quite close to the good fit obtained 
in the LO analysis nDS. Interestingly, the \ 2 °f the NLO fit of nDS is slightly smaller 
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The calculated scale evolution (solid lines) of the ratio F^/F^ compared with 
data [13J for several fixed values of x. The inner error bars are the statistical ones, 
ones stand for the statistical and systematic errors added in quadrature. 
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than the LO ones, lending further support to the validity of the global analysis. 



Set 


Ref. 


Q 2 /GeV 2 


-Ndata 


-Aparams 


x 2 


X 2 /N 


X 2 /d.o.f. 


This work 




1.69 


514 


16 


410.15 


0.798 


0.824 


EKS98 


[5j 


2.25 


479 




387.39 


0.809 




HKM 


[7\ 


1.0 


309 


9 


546.6 


1.769 


1.822 


HKN 


[8\ 


1.0 


951 


9 


1489.8 


1.567 


1.582 


nDS, LO 


[9J 


0.4 


420 


27 


316.35 


0.753 


0.806 


nDS, NLO 


[9J 


0.4 


420 


27 


300.15 


0.715 


0.764 



Table 3: The goodness of the fits obtained in different global analyses. 



To demonstrate the remaining uncertainties in the nPDFs, we show in Fig. [10] the 
comparison between this work (solid), EKS98 \6\ (dashed), HKM [7] (dotted), HKN 
(long-dashed) and nDS (NLO) [9J (dot-dashed) sets. The ratios R v , R£, R% and 
Rp 2 are plotted for A = 40 at scales Q 2 = 2.25 and 100 GeV 2 . We choose Calcium 
here as there are both small-x and larger-x DIS data and DY data available for this 
nucleus. The lower one of the scales considered is the initial scale in the EKS98 set. 

As can be seen in Fig. (TTJ the quantitative main difference between the present 
analysis and EKS98 lies in the small-x behaviour of sea quark and gluon modifications. 
For the sea quarks, the difference is merely due to the different form of the fit functions 
chosen: in the present work, shadowing in R$, and thus also that in Rp 2 , levels off 
faster. Like in the original EKS98 framework, the very-small-x behaviour of R G at is 
tied to that of Rp (but indirectly, through restricting the limits of the free parameters 
controlling the antishadowing maximum), thus also the gluon shadowing saturates now 
faster than in EKS98, and hence we have also somewhat less antishadowing in gluons. 
Recall also the small difference in the initial scales here and in EKS98. In the region 
x ~ 0.02 — 0.03, where the ratios Rq are indirectly constrained by the NMC data in 
Fig. [9] the results from the present work and EKS98 are very similar. 

Regarding all sets, we first notice that in the mid/large-x region x>0.1 the ratios 
Rp 2 are almost identical, thanks to the constraints given by the DIS data for the x, 
Q 2 and A dependence of Rp 2 . Since in the large-x region, x>0.3 or so, valence quarks 
dominate Rp , also the ratios Ry from different sets agree nicely there. 

The role of the DY data in pinning down both Ry and Rjj in the small/mid- x 
region 0.01^ x ^0.3 can be concretely seen in the figure. In the HKM |7J analysis 
(dotted lines), the DY data was not included. As a result, the HKM fit suggested 
Rg > 1 at a; > 0.1, which in turn compensated the smallness of Ry(x ~ 0.1) (see 
the left panel) in reproducing Rp 2 . The main improvement from HKM to HKN [S] 
was the inclusion of the DY data in the fit. This translates into better constrains and 
a better agreement with EKS98 for Ry over the whole x-region and also for Rg at 
0.01 5 x ^ 0.1. The fact that the ratios Ry from different global analyses agree so nicely 
is quite reassuring, as it demonstrates that the average valence quark modifications can 



21 




'l0~ 4 10" 3 10" 2 10"' 10" 4 10" 3 10" 2 10"' 10" 4 10" 3 10" 2 10"' 10" 4 10" 3 10" 2 10"' 1 



X 

Figure 10: (Colour online) Comparison of different nPDF modifications for Ca: This work, 
EKS98 [5J[6], HKM [7j, HKN [8j and nDS (NLO) [9] are plotted at scales Q 2 = 2.25 and 100 
GeV 2 . Calcium is used here as it is fairly well constrained by the data. 

be pinned down in a manner which does not depend much on the specific form chosen 
for the fit functions. 

At x > 0.2, where valence quarks start to dominate the quark sector, sea quarks are 
not sufficiently constrained by either DIS or DY data - hence the large variations in 
Rg from set to set. This is the case also in the very-small- x region x<0.01, in the 
absence of sufficient data constraints there. Thus, the very-small-x behaviour of Rg is 
specific to the form of the fit function chosen. 

As can be seen in Fig. [TO], the nuclear gluon distributions in general are still quite 
badly constrained, resulting in large differences between the different sets. In the 
absence of data which would sufficiently stringently constrain the gluon modifications 
over a wide enough x-range, the results from the global fits are bound to depend on the 
form of the fit functions chosen. To demonstrate this, we replot the ratio yp^F 2 Sn / j^F^ 
in Fig. [TT] for the six smallest-x panels of Fig. [HI As can be seen here, the log Q 2 slopes 
of F 2 Sn /F 2 become flatter in HKN and HKM than those in the present analysis, EKS98 
and nDS. The reason for this can be seen from the ratios Rq at x £ 0.02 — 0.04 in Fig.fTOl 
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and from Eq. fj 1 ip : the larger Rq is relative to Rp , the faster is the Q 2 dependence 
of Rp . However, as commented in the previous section, the small-x NMC data which 
would give at least some constraints for the gluons at x ~ 0.02 — 0.04, has a relatively 
small weight in the global analysis. All this makes it difficult to pin down the nuclear 
gluon modifications. 
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Figure 11: (Colour online) Comparison of the results from this analysis (solid), 
EKS98 (dashed), HKM (dotted), HKN (long-dashed) and nDS (dot-dashed) for the ratio 
m F 2 n /j2 F 2- As in Fi g- El (6 first panels there), the data is from NMC p]. 



5 Error Analysis 

Next, to quantify the above discussion on the uncertainties, we proceed to the error 
analysis, one of the goals in the present paper. We do this by using the Hessian 
method, which is one of the standard methods in multiparameter analyses as it takes 
the parameter correlation into account. The error matrix, or the Hessian matrix, is 
the inverse of the second derivative matrix of the fitting function \ 2 with respect to its 
free parameters. The Minuit fitting routine provides also this matrix along with the fit 
parameters [37J . Denoting the set of fit parameters by £ and the Hessian error matrix 
by H, the fitting function \ 2 can be expanded around the minimum £ as (See e.g. Ref. 
[38] . here we follow the notation of Ref. [39] ) 

Ax 2 = X\l + SO - x\l) = E H iMMr (12) 

i,3 
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The uncertainty of the fitted function F(x,£) is then 



[6F(x,i)] 2 = A x 2 J2 



) 




(13) 



assuming linear error propagation. However, the confidence region of a multivariable fit 
is different than that of a single variable fit and needs to be evaluated. The confidence 
level P of the normal distribution with N degrees of freedom can be written as 



where T(n) is the Gamma function. For one-parameter fit the one-cx error range results 
confidence level P = 0.6826 and Ax 2 = 1. Requiring the same confidence level for iV 
parameters one can now calculate the A% 2 . For example, for N = 16 one obtains 
A X 2 = 18.11. 

The error limits obtained using this method for the fit with the 16 free parameters 
in Table [2] are shown by the dashed lines in Fig. [12] for the ratios Ry, Rg, Rq and 
Rp 2 in the case of a Lead nucleus, A = 208. As can be seen from Fig. [121 an d as 
expected on the basis of Sec. HJ the ratio Ry is relatively well constrained. Also Rp 2 is 
rather well under control. At large x its errors naturally follow the small errors of Ry, 
thanks to the DIS data available. In the small/mid-x region both the DIS data and 
the DY data are necessary to pin down Ry and Rg. Towards smaller values of x the 
errors in Rg get larger due to the lack of high-precision constraints for the sea quarks 
there, but are nevertheless still constrained. As discussed in Sec. HJ the small-x errors 
of Rg shown, and thereby those in R§ , are specific to the small- a; behaviour assumed. 
Hence, the error bars given here are to be considered as lower limits. 

For the gluons, the very-small-x errors become quite large as there are no data 
constraints there to guide us. Similarly to the sea quark case, the error bars on gluon 
shadowing are fit function specific, and hence lower limits. However, as noticed in 
EKS98 and originally in Ref. [10] gluons do get somewhat better constrained at x ~ 
0.02 — 0.04, thanks to the NMC data. Note that the zero-error we obtain at the peak 
of the gluon antishadowing bump is an artifact due to the interplay between the free 
parameters and the momentum sum rule. 

To get physically more relevant estimates on the sea quark and gluon uncertainties 
for the mid- and large-x regions, we do the following. We free the parameters y a , p Va , 
Uei Py e and P (which control the magnitudes of the modifications in Rg and Rq) while 
keeping the location parameters x a and x e as well as the parameters controlling the 
small-x behaviour fixed to the values quoted in Table [2J Minimization of x 2 first with 
the freed sea quark parameters, then with the freed gluon parameters results in the 
wide bands shown by the dotted lines in Fig. [121 This demonstrates clearly how badly 
the nuclear sea quark and gluon modifications are constrained in the large-x region. 
Similar results have been presented before by the HKN group. Thus, as the error 




(14) 



24 



A = 208, Ql = 1.69 GeV 2 
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Figure 12: (Colour online) Fit errors at the initial scale Qq = 1.69 GeV 2 for Lead, shown by 
the dashed lines. For large- x sea and gluon modifications the errors shown by the dotted lines 
were calculated separately, see the text. The shaded (yellow on-line) band is the total error 
estimate obtain, see the text. The corresponding EKS98 results, evolved downwards from 
Qoeks = 2-25 GeV 2 , are shown by the dot-dashed (red) lines. An example of a stronger 
gluon shadowing is shown by dense-dashed (green) line. 

estimates for the present analysis, we give the shaded (yellow on-line) bands of the 
small-x and large-x errors, denoting them by "total errors" in Fig. [12j 

In Fig. [j~2"1we also show the comparison with the EKS98 modifications, evolved from 
a higher initial scale, Qqeks = 2-25 GeV 2 , down to the present one, Qq = 1.69 GeV 2 . 
Within the errors estimated, we can safely conclude that the old EKS98 parametriza- 
tion is fully consistent with the present ^-minimization analysis. As discussed in the 
previous section, the fact that EKS98 sea quarks and gluons lie somewhat below the 
results from this work, is mainly due to the different functional forms assumed for the 
fit functions at small values of x. We thus conclude that there is no need for releasing 
a new LO parametrization, since EKS98 still works very well. 
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6 Stronger gluon shadowing? 



Similarly to our earlier work EKS98, the present analysis suggests that the nuclear 
gluon modifications in the region x ~ 0.02 — 0.04 should be rather small, while the 
amounts of shadowing and thus antishadowing are much more weakly constrained. As 
the final task in this paper we discuss the possibility of a stronger gluon shadowing. 
Our main motivation for doing this is the inclusive charged-hadron data taken from 
D+Au collisions at RHIC by the BRAHMS collaboration [4T] , and the computation of 
the corresponding p? spectra in Ref. [42] using the strong gluon shadowing suggested in 
Refs. [21| 143] [12]. These data are advocated as a hint that a parton saturation regime 
could have been reached at RHIC [44J, so the degree of agreement with a DGLAP 
approach is of special interest. 

We construct our strong gluon shadowing example by changing only the parameter 
y a for the Carbon reference nucleus in Rq. Then, as seen in Fig. [12] the changes in the 
region x ~ 0.02 — 0.04 remain small but the amounts of antishadowing and (through 
momentum conservation) shadowing change. Increasing y a from 1.071 to 1.2 deepens 
the saturation level of gluon shadowing in Lead considerably, from 0.7 to 0.26. At the 
same time, the goodness x 2 /N of the overall fit weakens only slightly, from 0.80 to 
0.95, even if no x 2 minimization was performed. 

With the gluon shadowing much stronger than that of sea quarks, the log Q 2 slopes 
of Rp 2 at small x are initially negative. At the same time, due to the stronger gluon 
antishadowing, the scale evolution of Rg near x ~ 0.1 is slightly speeded up. These 
effects can be verified in Fig. [13] (compare with Fig. [2]). In fact, the latter effect is 
responsible for the deterioration of the goodness. We stress, however, that for this 
strong gluon shadowing example we have kept the quark sector as given in Table [2] 
After minimization, the changes in x 2 /N would become even smaller, demonstrating 
the fact that quite large changes in the gluon sector induce only small changes in the 
global x 2 • This is interesting when compared with the results of de Sassot and Florian 
j9], who get considerably worse x 2 values for stronger gluon shadowing. Apparently, 
the form of their fit is such that stronger gluon shadowing in small-x affects in the 
region x ~ 0.01 — 0.1 as well, thus changing the fit there. 

In Fig. [J4] we show the ratio Rdau for minimum bias single hadron production, 
defined as 

1 da DAu 

E> _ A dp T dr] /-, p-n 

Rdau - daPP , (15 ) 

dpxdri 

where pr and 77 are the hadronic transverse momentum and pseudorapidity, corre- 
spondingly. The BRAHMS data in the top panels are for Rr>Au(h + + h~) and in the 
bottom panels for i? D A u (^~)- The generic structure of the lowest order pQCD cross 
sections is given by 

a AB^ h+ x = ^2 ff( Xlj Q) ® fffa q) 3 a i+i^+i g, Dk ^ h+x (z,Qfl (16) 

ijkl 
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Figure 13: Scale evolution of the ratios Rg , Rq and Rp 2 for Carbon and Lead in the case of 
the strong gluon shadowing example considered in Fig. [T2j Notice the initial negative log Q 2 
slopes of Rg and hence also Rp 2 at small values of x. 
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where h is the hadron type, k labels the parton type, AB = DAu, pp and Dk->h+x(z, Qf) 
are the fragmentation functions at a fractional energy z = Eh/Ek and a factorization 
scale Qf. Detailed formulation of the computation can be found e.g. in [15] . Here 
we choose Q as the transverse momentum of the parton and Qf as the transverse 
momentum of the hadron. We use the KKP fragmentation functions [46] and the 
CTEQ6L1 free proton PDFs. We do not make attempt to correct for the fact that 
the KKP fragmentation functions correspond to the average h + + h~, even though the 
forward-rapidity data is for negative hadrons only. 
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Figure 14: (Colour online) Minimum bias inclusive hadron production cross sections in d+Au 
collisions divided by that in p+p collisions at \fs NN = 200 GeV at RHIC. The ratio -RdAu is 
shown as a function of hadrons transverse momentum at four different pseudorapidities. The 
BRAHMS data [41] are shown with the statistical error bars and the shaded systematic error 
limits. A pQCD calculation for h + + h~ production with the EKS98 nuclear modifications 
and KKP fragmentation functions is shown by the solid lines (red) and that with the strong 
gluon shadowing by the dashed lines (green). 

At small pseudorapidities, where both quark and gluon-initiated processes are im- 
portant, the stronger gluon antishadowing induces only a small correction to -Rdau 
but in a manner that the overall shape of the computed -Rdau agrees better with the 
BRAHMS data. At large pseudorapidities, corresponding to smaller x 2 , gluons be- 
come dominant. As discussed in [45], hadron production at, say, 1.5 GeV is biased to 
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partons at px ~ 3 GeV. Since %i = ^=(e n + e 2/2 ), small values of x 2 of the order 
0.001, start to play a role at 77 = 3. Integration over (or 22) however, smears the 
effects of the nuclear modifications which is why we do not see a larger change in -Rdau 
with the stronger gluon shadowing example considered. As shown in Ref. [42J, even 
more dramatic small- 2 behaviour of gluons, such as suggested in J2Tj ffi2 [12] , would 
obviously be needed to account for the BRAHMS data. Whether gluons with such 
shadowing, supplemented perhaps with stronger shadowing for the sea quarks as well, 
would maintain the good global fit to the DIS and DY data now obtained, remains to 
be seen. At the same time, dependence of the fragmentation functions on the hadron 
charge (negatives instead of the average h + + h~), should be studied in more detail 
within a consistent DGLAP framework. 

Due to the double integrations in computing the cross sections in Eq. [TBJ inclusion 
of the RHIC data for -Rdau in the global analysis is beyond the scope of the present 
paper. As further data constraints are absolutely necessary for pinning down the 
nuclear gluons, these data, in spite of their relatively large systematic errors, motivate 
us to do this in future. 

7 Summary 

In this study we have performed a global leading-order DGLAP analysis of the nPDFs 
using the EKS98 framework introduced in [5j [6]. Motivated by our previous work, 
we have introduced a piece-wize parametrization for the nuclear effects in the PDFs. 
Originally, the fit functions contained altogether 42 parameters. With the help of mo- 
mentum and baryon number conservation and the experience from EKS98, we reduce 
the number of relevant fit parameters down to 16. A best fit to the nuclear DIS and 
DY data was searched for this set of parameters through automated minimization of x 2 
using the Minuit program [37]. As a result, a very good fit to the iV = 514 data points 
at Q 2 > 1.69 GeV 2 was found, giving x 2 /N = 0.789 (or x 2 /d.o.f. = 0.82). No essential 
improvement over EKS98 was found, however, as the EKS98 modifications lead to an 
equally good fit quality, x 2 /N = 0.809 (for N = 479 datapoints at Q 2 > Qo,eks9s)- 

Relative to the old EKS98, the present analysis suggests slightly less shadowing for 
the gluons and sea quarks. This, however, is merely due to the different forms of the 
fit functions adopted in the region where no stringent constraints from the data are 
available. We also compared the obtained nuclear effects to those obtained by other 
global analyses, HKM, HKN, and nDS. The valence quark modifications do not deviate 
much from one set to another but the smallest-x and large-x modifications of gluons 
and sea quarks differ in a major way. This reflects the fact that especially the nuclear 
gluons are badly constrained in these regions. 

To quantify the uncertainties in our analysis, we obtained the error estimates by 
using the Hessian method based on the information given by Minuit. The error esti- 
mates obtained also nicely further confirm the validity of EKS98, as it is shown to be 
fully consistent with the present analysis. 
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To get a hold on the uncertainties in the large- x regions of gluons and sea quarks, 
we computed the large-x errors separately. These, considered together with the small- x 
errors on the best fit confirm the conclusions from the comparison between different 
analyses: the valence quark distributions are relatively well, and independently from 
the fit-function form, constrained over the whole x region. For the sea quarks, the large- 
re (x^0.3) errors become very large, and for the small-x behaviour clearly depends 
on the fit function form. For gluons, our analysis shows that presently one can to 
some extent constrain the gluons in the region x ~ 0.02 — 0.04 but hardly at all 
in the large-x region, and only in a fit-function-dependent manner at small x through 
momentum conservation. We also note that the relatively small error estimate obtained 
at x ~ 0.02 — 0.04 for gluons may depend somewhat on the framework chosen, as the 
gluon fit parameters were drifting to the limits imposed. This obviously leaves room 
for further improvements in the future. An obvious further improvement of the present 
analysis is its extension to NLO. 

As the DIS and DY data are not able to stringently pin down the gluon modifica- 
tions, further constraints are obviously needed. In thinking of possible additional data 
sets to be included in the global analysis in the future, we considered an example of 
a stronger gluon shadowing without doing a x 2 minimization. First, we showed that 
quite large variations in the gluon modifications can be absorbed in the quark sector 
and thus hidden by the good \ 2 values obtained. Then, motivated by Ref. [42], we 
computed the nuclear modification ratio -Rdau of inclusive hadron production in d+Au 
relative to that in pp, using both the EKS98 modifications and the strong gluon shad- 
owing example. Comparisons against the BRAHMS data j4T] here and in Ref. [42] lend 
support to more shadowed gluons than in the present EKS98 framework. At RHIC, the 
d+Au data is evidently very valuable for getting further constraints for nuclear gluons 
in particular. This in turn demonstrates the importance of running a parallel p+Pb 
program at the LHC, where pQCD factorization and nPDFs could be tested further in 
a wide range of x and Q 2 . 
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